Geophysical Acquisition Design
Three-dimensional (3D) seismic surveys in the exploration industry can cost tens to hundreds of million dollars. 3D survey data are expensive to acquire because of the necessary short receiver and source intervals, which are primarily driven by signal and noise issues, and the wide range of offsets and azimuths that are required by imaging methods. Ultimately the quality of the processed data and subsequent interpretations are dependent on how well the survey was designed and acquired. No matter how sophisticated the methods used to process the data are, the processing cannot compensate for a poorly planned 3D survey and inadequate sampling, which can end up making the survey a waste of time and resources if it does not meet its objectives.
Acquisition design typically starts by assuming different designs and then testing and evaluating each design. The first step is a survey geometry analysis based strictly on surface bins. For each survey geometry, bin-based attributes, such as the fold and the offset and azimuth distributions, are computed assuming uniform earth properties. The attributes are used to select a few survey geometries for further consideration. Then if the geometry is complex, a geological model may be constructed and used for further evaluation of a few designs selected based on the geometry analysis. Source and receiver locations are assumed and seismic waves are then propagated from each of the sources to the target and then to each of the receivers using some form of computer simulation. The computer simulation may be done with ray-tracing methods or by solving the one-way or two-way wave-equation. Wave-equation methods are more accurate for complex media but are much more computer intensive. Methods differ in how they use and evaluate the computer simulations for each survey design.
“Brute-Force” Comparison Methods
Comparison of two survey designs using a 3D geological model can be performed by fully simulating, and then, separately processing and imaging the computer simulated data corresponding to the two designs. This may be called a “brute-force” method. The two imaged volumes are compared to determine whether or not they meet the survey objectives. This typically requires wave-equation simulation with finite-difference methods to generate data that most accurately approximate the data that would be recorded in the field. Then each dataset must be separately processed and imaged as would be the field data, typically with finite-difference migration methods such as Reverse Time Migration (RTM). The wave equation simulation and processing are both time consuming and expensive; it can be prohibitively expense to carry out this effort for several models and for many different acquisition geometries. Most commonly this step is done after other analysis to compare only two scenarios.
There are a number of limitations to the brute-force methods (full simulation and comparison) in addition to the cost and time required. Developing a metric for quantifying the difference between the two imaged volumes from different design scenarios is difficult and often comes down to a “beauty contest.” An interpreter may look at the images and pick one as being more “geological”. Seismic attributes such as an amplitude map at a subsurface layer, may be computed and compared. The choice of the best geometry is subjective because often one area looks better with one design and another area looks better for the other design. A second problem is that the geological model is uncertain, and there is no way to validate the sensitivity to changes in the model. It is prohibitively expensive to compare two full 3D surveys for two or more models. Typically, models are built from existing information such as processed spec data and known geological information in the area. Often there are several models in an area, particularly in exploration settings, suggesting many possible geological scenarios and plays. A third problem with this approach is that this and other approaches do not interrogate the crucial and fundamental need to obtain an accurate velocity model for the imaging. Typically, a single smoothed version of the known model velocities is used for imaging both of the two designs. This ignores the fact that the imaging, and in particular the imaged amplitudes, critically depend on the velocities used for imaging, and these velocities are unknown and must be determined from the data. The method does not evaluate whether or not one survey design is better than another survey design for velocity determination. Conventional velocity determination is a tedious process that requires subjective input and parameterization, and thus, it is difficult to compare the results of two designs because there can be differences from having two different people doing the determination or from one person doing it twice having learned from the first set.
The brute-force method (full simulation and comparison) involves shooting two geometries with the computer instead of in the field, processing the two geometries as would be done with field data, and then comparing the outputs. One example is illustrated in the publications of Carl Regone (28) comparing narrow and wide azimuth survey geometries for subsalt imaging. Another is the publication of L'Heureux (29) comparing receivers in a marine streamer to those on the ocean bottom. A third example is disclosed by Rekdal (34) who evaluates the effect of holes in the data coverage of a marine survey by observing the data quality degradation for processed and imaged data of a simulated survey with holes compared to the full-coverage survey, in which “data quality degradation” is related to the observation of artifacts in the images. Rekdal assumes a velocity model for imaging and does not solve for a velocity model or use a quantitative measure related to the model for evaluation.
A fourth example does test the two geometries for velocity determination, but only in 2D because of the cost and time required. Baptiste et al., [12] use 2D full waveform inversion to investigate survey parameters for an ocean bottom cable acquisition problem with increasing source-to-receiver offsets. Although their approach uses full waveform inversion in survey design, they are performing it in the same manner as one would do if processing the field data, and each geometry is processed and inverted separately. They restrict the test to 2D because of speed issues.
Attribute Methods
Because of the cost and time associated with simulating and processing 3D data, it is more common to compute and display attributes related to illumination qualities of the acquisition geometry instead of fully processing and imaging the simulated data. This may be called the “attribute method.” The attributes may be computed for a reflective interface in the model or computed for a single target location. Typical attributes may include illumination hits per grid cell, simulated migration amplitude maps, and angle coverage displays. The simulation method may be a simple ray-tracing method or a more sophisticated wave-equation method. In the latter case, it is more common to compute only attributes for one target location because of the increase in cost. The problem with all the attribute methods is that it is not clear how the attributes are linked to the survey objectives and the ability to image a specific area of the subsurface. It can be difficult to evaluate cost/benefit tradeoffs between different geometries for the purposes of reservoir delineation and characterization. It is typical to display an image of the attributes and pick the “better” pattern without any direct relationship between the pattern and the objectives. In addition, attributes related to the imaged amplitude or energy are used, but for seismic waves, amplitudes are secondary to kinematics based on the traveltimes of the seismic waves. Imaging requires good imaging velocities determined from the traveltimes. None of these methods investigate the ability to determine the velocities for the imaging in the first place. Some examples of the attribute approach are briefly described in the following sections.
A commercial package for performing ray tracing and attribute analysis is provided by Norsar (www.norsar.com). The package traces rays for the geometrical model and then computes illumination attributes and imaged amplitudes.
Campbell [30] teaches a method to design a targeted survey based on an illumination analysis in which emergent ray paths are observed, and a determination is made to put sources and receivers where the rays are concentrated.
Bernitas [31] teaches a hybrid wave-equation ray-tracing method for survey design that computes and analyzes wave function amplitudes and illumination images.
Velduizen[32] uses a focal beam approach to predict the imaged amplitudes at the target. For a single target location, she determines the range of illumination angles and azimuths for which the amplitude is uniform.
Cao [33] teaches a use of the two-way wave equation and computing and analyzing the energy and associated energy attributes for use in survey design.
Optimization Methods
Instead of comparing attribute displays, a number of papers propose to invert for the best geometry that minimizes an objective function based on an attribute. This may be called the “optimization method.” Mauerer et al. [1] provide a comprehensive review of several survey design methods that are posed as inverse problems. The goal is to determine optimal survey parameters such as source and receiver locations for geophysical acquisition to address specific objectives of “better” illuminating subsurface models. References that pose the survey design as an optimization problem include the following.
Mauerer and Boerner [2] pose the survey design as a statistical inverse problem and use nonlinear optimization method to solve for survey parameters.
Curtis et al., [4, 5, 6, 7], solve a borehole acquisition design problem as a Bayesian survey design problem as well as maximizing the eigenvalues of inverse operators. Both approaches are dependent on posing it as inverse problem. Haber et al. [9], Coles and Morgan [3] also apply a Bayesian approach to solve survey design problems.
Routh et al. [10] and Oldenborger and Routh [11] maximize the resolution measure of a point spread function to obtain optimal survey parameters by solving a highly nonlinear optimization problem for an electromagnetic and electrical problem. Ajo-Franklin [8] extended the Routh at al. [10] approach to consider many point spread functions to solve for survey design parameters for a borehole problem. All of the approaches are formulated as a nonlinear inverse problem.
All of these optimization methods pose the survey design as an optimization problem, and it is well recognized that such approaches are extremely nonlinear and difficult to solve. Moreover, such optimization methods can rarely be applied to large scale experiments such as full scale 3D seismic surveys. In the concluding remarks of their review paper, Mauerer et al. [1] state that most survey design algorithms are computationally expensive, which imposes severe constraints on applications to realistic design problems. As a future initiative, they suggest improving the computational efficiency of key modules such as forward simulation engines that are called many times while solving a design algorithm. Even with such advances, the method is unlikely to be fast or robust enough to be used for realistic large scale 3D survey design problems commonly encountered in seismic industry. In the approach described in the present disclosure, the whole premise of posing the survey design as an optimization problem is avoided. That not only avoids solving a complicated nonlinear optimization problem but also allows one to consider realistic large scale survey design problems commonly encountered in seismic industry using a simultaneous source method.
Full Wavefield Inversion with Simultaneous Sources
Full-wavefield inversion of surface seismic data is an emerging technology for processing seismic surveys; however, even with modern computing power, seismic full wavefield inversion is still a computationally expensive endeavor. However, the benefit of obtaining a detailed representation of the subsurface using this method provides strong incentive to overcome this impediment. Development of algorithms and workflows that lead to faster turn around time is a key step towards making this technology feasible for field scale data. Seismic full waveform inversion involves several iterations of forward and adjoint simulation of the data. Therefore techniques that reduce the cost of forward and adjoint computation runs will allow users to solve larger-scale problems in a reasonable amount of time.
Geophysical inversion [13, 14] attempts to find a model of subsurface properties that optimally explains observed data and satisfies geological and geophysical constraints. There are a large number of well-known methods of geophysical inversion. These well-known methods fall into one of two categories, iterative inversion and non-iterative inversion. The following are definitions of what is commonly meant by each of the two categories:
Non-iterative inversion—inversion that is accomplished by assuming some simple background model and updating the model based on the input data. This method does not use the updated model as input to another step of inversion. For the case of seismic data these methods are commonly referred to as imaging, migration, diffraction tomography or Born inversion.
Iterative inversion—inversion involving repetitious improvement of the subsurface properties model such that a model is found that satisfactorily explains the observed data. If the inversion converges, then the final model will better explain the observed data and will more closely approximate the actual subsurface properties. Iterative inversion usually produces a more accurate model than non-iterative inversion, but is much more expensive to compute.
Iterative inversion is generally preferred over non-iterative inversion, because it yields more accurate subsurface parameter models. Unfortunately, iterative inversion is so computationally expensive that it is impractical to apply it to many problems of interest. This high computational expense is the result of the fact that all inversion techniques require many computationally intensive simulations. The compute time of any individual simulation is proportional to the number of sources to be inverted, and typically there are large numbers of sources in geophysical data, where the term source as used in the preceding refers to an activation location of a source apparatus. The problem is exacerbated for iterative inversion, because the number of simulations that must be computed is proportional to the number of iterations in the inversion, and the number of iterations required is typically on the order of hundreds to thousands.
The most commonly employed iterative inversion method employed in geophysics is cost function optimization. Cost function optimization involves iterative minimization or maximization of the value, with respect to the model M, of a cost function S(M) which is a measure of the misfit between the calculated and observed data (this is also sometimes referred to as the objective function), where the calculated data are simulated with a computer using the current geophysical properties model and the physics governing propagation of the source signal in a medium represented by a given geophysical properties model. The simulation computations may be done by any of several numerical methods including but not limited to finite difference, finite element or ray tracing. The simulation computations can be performed in either the frequency or time domain.
Cost function optimization methods are either local or global [15]. Global methods simply involve computing the cost function S(M) for a population of models {M1, M2, M3, . . . } and selecting a set of one or more models from that population that approximately minimize S(M). If further improvement is desired this new selected set of models can then be used as a basis to generate a new population of models that can be again tested relative to the cost function S(M). For global methods each model in the test population can be considered to be an iteration, or at a higher level each set of populations tested can be considered an iteration. Well known global inversion methods include Monte Carlo, simulated annealing, genetic and evolution algorithms.
Unfortunately global optimization methods typically converge extremely slowly and therefore most geophysical inversions are based on local cost function optimization. Algorithm 1 summarizes local cost function optimization.                1. selecting a starting model,        2. computing the gradient of the cost function S(M) with respect to the parameters that describe the model,        3. searching for an updated model that is a perturbation of the starting model in the negative gradient direction that better explains the observed data.        
Algorithm 1—Algorithm for Performing Local Cost Function Optimization.
This procedure is iterated by using the new updated model as the starting model for another gradient search. The process continues until an updated model is found that satisfactorily explains the observed data. Commonly used local cost function inversion methods include gradient search, conjugate gradients and Newton's method.
Local cost function optimization of seismic data in the acoustic approximation is a common geophysical inversion task, and is generally illustrative of other types of geophysical inversion. When inverting seismic data in the acoustic approximation the cost function can be written as:
                              S          ⁡                      (            M            )                          =                              ∑                          g              =              1                                      N              g                                ⁢                                          ⁢                                    ∑                              r                =                1                                            N                r                                      ⁢                                                  ⁢                                          ∑                                  t                  =                  1                                                  N                  t                                            ⁢                                                          ⁢                              W                ⁡                                  (                                                                                    ψ                        calc                                            ⁡                                              (                                                  M                          ,                          r                          ,                          t                          ,                                                      w                            g                                                                          )                                                              -                                                                  ψ                        obs                                            ⁡                                              (                                                  r                          ,                          t                          ,                                                      w                            g                                                                          )                                                                              )                                                                                        (                  Eqn          .                                          ⁢          1                )            where:    S=cost function,    M=vector of N parameters, (m1, m2, . . . mN) describing the subsurface model,    g=gather index,    wg=source function for gather g which is a function of spatial coordinates and time, for a point source this is a delta function of the spatial coordinates,    Ng=number of gathers,    r=receiver index within gather,    Nr=number of receivers in a gather,    t=time sample index within a trace,    Nt=number of time samples,    W=minimization criteria function (a common choice is W(x)=x2, which is the least squares (L2) criteria),    ψcalc=calculated seismic pressure data from the model M, often expressed as a matrix,    ψobs=measured seismic pressure data.
The gathers can be any type of gather that can be simulated in one run of a seismic forward modeling program. Usually the gathers correspond to a seismic shot, although the shots can be more general than point sources. For point sources the gather index g corresponds to the location of individual point sources. For plane wave sources g would correspond to different plane wave propagation directions. This generalized source data, ψobs, can either be acquired in the field or can be synthesized from data acquired using point sources. The calculated data ψcalc on the other hand can usually be computed directly by using a generalized source function when forward modeling. For many types of forward modeling, including finite difference modeling, the computation time needed for a generalized source is roughly equal to the computation time needed for a point source.
Equation (1) can be simplified to:
                              S          ⁡                      (            M            )                          =                              ∑                          g              =              1                                      N              g                                ⁢                                          ⁢                                    W              ⁡                              (                                  δ                  ⁡                                      (                                          M                      ,                                              w                        g                                                              )                                                  )                                      .                                              (                  Eqn          .                                          ⁢          2                )            where the sum over receivers and time samples is now implied and,δ(M,wg)=ψcalc(M,wg)−ψobs(wg)  (Eqn. 3)
Inversion attempts to update the model M such that S(M) is a minimum. This can be accomplished by local cost function optimization which updates the given model M(k) as follows:M(k+1)=M(k)−α(k)∇MS(M)  (Eqn. 4)
where k is the iteration number, a is the scalar size of the model update, and ∇MS(M) is the gradient of the misfit function, taken with respect to the model parameters. The model perturbations, or the values by which the model is updated, are calculated by multiplication of the gradient of the objective function with a step length a, which must be repeatedly calculated.
From equation (2), the following equation can be derived for the gradient of the cost function:
                                          ∇            M                    ⁢                      S            ⁡                          (              M              )                                      =                              ∑                          g              =              1                                      N              g                                ⁢                                          ⁢                                                    ∇                M                            ⁢                              W                ⁡                                  (                                      δ                    ⁡                                          (                                              M                        ,                                                  w                          g                                                                    )                                                        )                                                      .                                              (                  Eqn          .                                          ⁢          5                )            
So, to compute the gradient of the cost function, one must separately compute the gradient of each gather's contribution to the cost function, then sum those contributions. Therefore, the computational effort required for computing ∇MS(M) is Ng times the compute effort required to determine the contribution of a single gather to the gradient. For geophysical problems, Ng usually corresponds to the number of geophysical sources and is on the order of 10,000 to 100,000, greatly magnifying the cost of computing ∇MS(M).
Note that computation of ∇MW(δ) requires computation of the derivative of W(δ) with respect to each of the N model parameters mi. Since for geophysical problems N is usually very large (usually more than one million), this computation can be extremely time consuming if it had to be performed for each individual model parameter. Fortunately, the adjoint method can be used to efficiently perform this computation for all model parameters at once [13]. The adjoint method for the least squares objective function and a gridded model parameterization is summarized by the following algorithm:                1. Compute forward simulation of the data using the current model and the gather signature wg as the source to get ψcalc(M(k),wg),        2. Subtract the observed data from the simulated data giving δ(M(k),wg),        3. Compute the reverse simulation (i.e. backwards in time) using δ(M(k),wg) as the source producing ψadjoint(M(k),wg),        4. Compute the integral over time of the product of ψcalc(M(k),wg) and ψadjoint(M(k),wg) to get ∇MW(δ(M,wg))—where ψadjoint is the adjoint matrix for matrix ψcalc.        
Algorithm 2—Algorithm for Computing the Least-Squares Cost-Function Gradient of a Gridded Model Using the Adjoint Method.
While computation of the gradients using the adjoint method is efficient relative to other methods, it is still very costly. In particular, the adjoint method requires two simulations, one forward in time and one backward in time, and for geophysical problems these simulations are usually very compute intensive. Also, as discussed above, this adjoint method computation must be performed for each measured data gather individually, increasing the compute cost by a factor of Ng.
The compute cost of all categories of inversion can be reduced by inverting data from combinations of the sources, rather than inverting the sources individually. This may be called simultaneous source inversion. Several types of source combination are known including: coherently sum closely spaced sources to produce an effective source that produces a wavefront of some desired shape (e.g. a plane wave), sum widely spaces sources, or fully or partially stacking the data before inversion.
The compute cost reduction gained by inverting combined sources is at least partly offset by the fact that inversion of the combined data usually produces a less accurate inverted model. This loss in accuracy is due to the fact that information is lost when the individual sources are summed, and therefore the summed data does not constrain the inverted model as strongly as the unsummed data. This loss of information during summation can be minimized by encoding each shot record before summing. Encoding before combination preserves significantly more information in the simultaneous source data, and therefore better constrains the inversion [16, 17]. Encoding also allows combination of closely spaced sources, thus allowing more sources to be combined for a given computational region. Various encoding schemes can be used with this technique including time shift encoding and random phase encoding. The remainder of this Background section briefly reviews various published geophysical simultaneous source techniques, both encoded and non-encoded.
Van Manen [19] suggests using the seismic interferometry method to speed up forward simulation. Seismic interferometry works by placing (in the simulation) sources everywhere on the boundary of the region of interest. These sources are modeled individually, and the wavefield at all locations for which a Green's function is desired is recorded. The Green's function between any two recorded locations can then be computed by cross-correlating the traces acquired at the two recorded locations and summing over all the boundary sources. If the data to be inverted have a large number of sources and receivers that are within the region of interest (as opposed to having one or the other on the boundary), then this is a very efficient method for computing the desired Green's functions. However, for the seismic data case, it is rare that both the source and receiver for the data to be inverted are within the region of interest. Therefore, this improvement has very limited applicability to the seismic inversion problem.
Berkhout [20] and Zhang [21] suggest that inversion in general can be improved by inverting non-encoded simultaneous sources that are summed coherently to produce some desired wave front within some region of the subsurface. (The main purpose of encoding is to make the source summation incoherent [17].) For example, point source data could be summed with time shifts that are a linear function of the source location to produce a down-going plane wave at some particular angle with respect to the surface. This technique could be applied to all categories of inversion. A problem with this method is that coherent summation of the source gathers necessarily reduces the amount of information in the data. So for example, summation to produce a plane wave removes all the information in the seismic data related to travel time versus source-receiver offset. This information is critical for updating the slowly varying background velocity model, and therefore Berkhout's method is not well constrained. To overcome this problem, many different coherent sums of the data (e.g. many plane waves with different propagation directions) could be inverted, but then efficiency is lost since the cost of inversion is proportional to the number of different sums inverted. Herein, such coherently summed sources are called generalized sources. Therefore, a generalized source can either be a point source or a sum of point sources that produces a wave front of some desired shape.
Van Riel [22] suggests inversion by non-encoded stacking or partial stacking (with respect to source-receiver offset) of the input seismic data, then defining a cost function with respect to this stacked data which will be optimized. Thus, this publication suggests improving cost function based inversion using non-encoded simultaneous sources. As was true of the Berkhout's [6] simultaneous source inversion method, the stacking suggested by this method reduces the amount of information in the data to be inverted and therefore the inversion is less well constrained than it would have been with the original data.
Mora [23] proposes inverting data that is the sum of widely spaced sources. Thus, this publication suggests improving the efficiency of inversion using non-encoded simultaneous source simulation. Summing widely spaced sources has the advantage of preserving much more information than the coherent sum proposed by Berkhout. However, summation of widely spaced sources implies that the aperture (model region inverted) that must be used in the inversion must be increased to accommodate all the widely spaced sources. Since the compute time is proportional to the area of this aperture, Mora's method does not produce as much efficiency gain as could be achieved if the summed sources were near each other.
Ober [24] suggests speeding up seismic migration, a special case of non-iterative inversion, by using simultaneous encoded sources. After testing various coding methods, Ober found that the resulting migrated images had significantly reduced signal-to-noise ratio due to the fact that broad band encoding functions are necessarily only approximately orthogonal. Thus, when summing more than 16 shots, the quality of the inversion was unsatisfactory. Since non-iterative inversion is not very costly to begin with, and since high signal-to-noise ratio inversion is desired, this technique is not widely practiced in the geophysical industry.
Ikelle [25] suggests a method for fast forward simulation by simultaneously simulating point sources that are activated (in the simulation) at varying time intervals. A method is also discussed for decoding these time-shifted simultaneous-source simulated data back into the separate simulations that would have been obtained from the individual point sources. These decoded data could then be used as part of any conventional inversion procedure. A problem with Ikelle's method is that the proposed decoding method will produce separated data having noise levels proportional to the difference between data from adjacent sources. This noise will become significant for subsurface models that are not laterally constant, for example from models containing dipping reflectors. Furthermore, this noise will grow in proportion to the number of simultaneous sources. Due to these difficulties, Ikelle's simultaneous source approach may result in unacceptable levels of noise if used in inverting a subsurface that is not laterally constant.
Source encoding proposed by Krebs et al. [17], which is incorporated herein by reference in all jurisdictions that allow it, is a very cost effective method to invert full wave field data. (The same approach of simultaneous inversion of an encoded gather will work for receivers, either via source-receiver reciprocity or by encoding the actual receiver locations in common-source gathers of data.) For fixed receivers, the forward and adjoint computations need to be performed only for a single effective source. Given the fact that hundreds of shots are recorded for typical 2D acquisition geometries, and thousands in the case of 3D surveys, computational savings from this method are quite significant. In practice, a fixed-receiver assumption is not strictly valid for most common field data acquisition geometries. In the case of marine streamer data, both sources and receivers move for every new shot. Even in surveys where the locations of receivers are fixed, the practice often is that not all receivers are “listening” to every shot, and the receivers that are listening can vary from shot-to-shot. This also violates the “fixed-receiver assumption.” In addition, due to logistical problems, it is difficult to record data close to the source, and this means that near-offset data are typically missing. This is true for both marine and land surveys. Both of these factors mean that for a simultaneous source gather, every receiver location will be missing data for some source shots. In summary, in simultaneous encoded-source inversion, for a given simultaneous encoded gather, data are required at all receiver locations for every shot, and this may be referred to as the fixed-receiver assumption of simultaneous encoded-source inversion.
Routh et al. [18] propose an adjustment to straightforward application of simultaneous encoded sources (and/or receivers) inversion that would enhance its performance when the fixed-receiver assumption is compromised. Routh et al. [18] use a cross-correlation objective function to achieve the application of simultaneous encoded-sources for non-fixed data. Therefore many survey acquisition designs that are based on non-fixed spread such as narrow azimuth, wide azimuth, rich azimuth, coil geometry, ocean bottom cables, ocean bottom nodes and vertical seismic profiling can be investigated using the method described by Routh et al. [18]. Other published methods for mitigating the effects of non-listening receivers in simultaneous encoded source inversion include U.S. Pat. Nos. 8,437,998, 8,775,143 and U.S. Patent Publication No. 2013/0238246.